// The core optimisations file
// contains files that allow program to interface with the extensions

#include "optimisations.h"
#include <windows.h>
#include <stdio.h>
#include "control.h"
#include <cstdlib> 
#include <ctime> 
#include <math.h>
#include <time.h>
#include "heap.h"
#include "fheap.h"
#include <mvicon.h>
#include <config.h>
#include "save.h"


// variables required for panel drawing originaly in pvengine.cpp
double maxclrbackup; // GLOBAL VARIABLE
HWND                    render_window ;
HPALETTE                hPalApp ;
int                     render_bitmap_depth = -1 ;
long                    render_bitmap_bpl ;
uchar                   *render_bitmap_surface ;
uchar                   *render_priority_surface ;
BitmapInfo              render_bitmap ;
BitmapInfo              bitmap_template ;
HINSTANCE               hInstance ;
unsigned                render_width ;
unsigned                render_height ;
int                     renderwin_xoffset ;
int                     renderwin_yoffset ;

//bool cellset = false;
HANDLE heapMutex;

//setings variables
//Opts x;
bool optimisation1 = false;//resolution subdivision
bool optimisation2 = false;//poisson distibution
bool optimisation3 = false;
bool optimisation4 = false;
bool optimisation5 = true;
bool vicon = false;

//tracing variables
bool tracing = false;
int area;
int pixelsPerSec = 0;
int lastPixelsPerSec = 0;
time_t lasttime = 0;
int shortest;
int longest;
bool newview = false;
bool refreshView = false;
int maxDelay=1;

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Display
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//true = normal display, false equals delay display
bool screenType = true;
void toggleScreen()
{
	screenType=!screenType;
	InvalidateRect (render_window, NULL, FALSE);
}
bool getScreen()
{
	return screenType;
}
void dropScreen()
{
	if (tracing)
	{
		WaitForSingleObject( heapMutex, INFINITE);  // no time-out interval
		save();
		ReleaseMutex(heapMutex);// release mutex
	}
}
bool working = false;
bool getWorking()
{
	return working;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Pixel representation and data structure
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// the struct will hold all the information about the pixels
struct Cell{
	Cell():trace(TraceType::none),area(1),delay(0){}
	TraceType trace;
	int area;
	int delay;
};
// the largest area possible
int maxarea;
//the number of pixels at area
int pixOffset[maxsize];
// a table that contains the info for each pixel
Cell pixelInfo [maxsize][maxsize];
Cell tempInfo [maxsize][maxsize];
void backupInfo()
{
	memcpy(tempInfo, pixelInfo, sizeof(tempInfo));
}
// a fibonacci heap that holds the queue of pixels to render
Heap *minQueue;
HeapD<FHeap> heapD;
bool heapSetup = false;
bool initialiseFibHeap()
{
	if (!heapSetup)
	{
		heapMutex = CreateMutex( 
			NULL,              // default security attributes
			FALSE,             // initially not owned
			NULL);             // unnamed mutex
		if (heapMutex == NULL) 
		{
			printf("CreateMutex error: %d\n", GetLastError());
			return false;
		}
		heapSetup = true;
	}
	//aquire mutex
	WaitForSingleObject( heapMutex, INFINITE);  // no time-out interval
	minQueue = heapD.newInstance(render_width*render_height);
	ReleaseMutex(heapMutex);// release mutex
	return heapSetup;
}


PixelRender getQueue()
{
	//aquire mutex
	WaitForSingleObject( heapMutex, INFINITE);  // no time-out interval
	PixelRender cell;
	if (minQueue->nItems()==0)
	{
		ReleaseMutex(heapMutex);
		return cell;
	}
	//get item from queue
	int id = minQueue->deleteMin();
	cell.x = id%render_width;
	cell.y = id/render_width;
	cell.trace = pixelInfo[cell.y][cell.x].trace;
	cell.area = pixelInfo[cell.y][cell.x].area;
	//reset cell
	pixelInfo[cell.y][cell.x].trace=TraceType::none;
	ReleaseMutex(heapMutex);
	return cell;
}
bool addToQueue(int x,int y,TraceType trace,int area,long delay,bool serialise)
{
	//check we are actually adding something
	if (trace==2)
		return false;
	//aquire mutex if requred
	if (serialise)
		WaitForSingleObject( heapMutex, INFINITE);  // no time-out interval
	//calculate value in queue

	priorityImage(x,y,area,delay);

	int i = (y*render_width) + x;
	//check if already in the queue
	if(pixelInfo[y][x].trace<=1)
	{
		//update pixel priority in queue
		if (pixelInfo[y][x].delay > delay)
		{
			pixelInfo[y][x].delay = delay;
			minQueue->insert(i,delay);
		}

		//if we are already set to trace then do not change to the less powerful silent trace
		if (!pixelInfo[y][x].trace==TraceType::trace)
			pixelInfo[y][x].trace=trace;

		//update quality if it is a larger area
		if (pixelInfo[y][x].area<area)
			pixelInfo[y][x].area;
	}
	else
	{
		//add to queue and set
		pixelInfo[y][x].trace=trace;
		pixelInfo[y][x].area=area;
		pixelInfo[y][x].delay=delay;
		minQueue->insert(i,delay);

	}
	//char out [20];
	//wsprintf(out, "\n%d", i);
	//OutputDebugString(out);
	if (serialise)
		ReleaseMutex(heapMutex);// release mutex
	return true;
}
void emptyQueue(bool emptyTable, bool serialise)
{
	//aquire mutex
	if (serialise)
		WaitForSingleObject( heapMutex, INFINITE);  // no time-out interval
	minQueue->~Heap();
	delete minQueue;
	minQueue = heapD.newInstance(render_width*render_height);
	for(int y = 0; y<render_height && emptyTable==true; y++)
	{
		for(int x = 0; x<render_width; x++)
		{
			//set to 0 if not optimisation 1
			pixelInfo [y][x].trace=TraceType::none;
			pixelInfo [y][x].area=0;
			pixelInfo [y][x].delay=0;
		}
	}

	if (serialise)
		ReleaseMutex(heapMutex);// release mutex
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											image comparison
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int maxPix;
bool similar [maxsize][maxsize];
int noSimilarPixels = 0;
int lastPercent = 0;
int imagesize=1;
uchar *copy_bitmap ;
bool compare = false;
bool allowCompare = false;
void resetCompare()
{
        if (!allowCompare)
                return;
        WaitForSingleObject( heapMutex, INFINITE);
        noSimilarPixels = 0;
        for (int y = 0; y < render_height; y++)
        {
                for (int x = 0; x < render_width; x++)
                {
                        uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x*3 ;
                        uchar* q = copy_bitmap + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x*3 ;
                        if ((*p++==*q++)&&(*p++==*q++)&&(*p++==*q++))
                        {
                                similar[y][x]=true;
                                noSimilarPixels++;
                        }
                        else
                        {
                                similar[y][x]=false;
                        }
                }
        }
        ReleaseMutex(heapMutex);
}
void updateRegion(int x, int y, int quality)
{
        if (!allowCompare||!compare)
                return;
        WaitForSingleObject( heapMutex, INFINITE);
        for (int y1 = y; y1 < y+quality && y1<render_height; y1++)
        {
                for (int x1=x; x1 < x+quality && x1<render_width; x1++)
                {
                        int offset = (render_bitmap.header.biHeight - 1 - y1) * render_bitmap_bpl + x1*3;
                        uchar* p = render_bitmap_surface + offset;
                        uchar* q = copy_bitmap + offset;
                        if (*p++==*q++&&*p++==*q++&&*p++==*q++)
                        {
                                if (!similar[y1][x1])
                                {
                                        similar[y1][x1]=true;
                                        noSimilarPixels++;
                                }
                        }
                        else
                        {
                                if (similar[y1][x1])
                                {
                                        similar[y1][x1]=false;
                                        noSimilarPixels--;
                                }
                        }
                }
        }
        if (lastPercent!=(int)(((double)noSimilarPixels*100)/maxPix))
                repaint(compare_value);
        ReleaseMutex(heapMutex);
}
void requiredSize(int x)
{
	if (x!=imagesize)
	{
		WaitForSingleObject( heapMutex, INFINITE);
		imagesize=x;
		noSimilarPixels=0;
		copy_bitmap = (BYTE *) calloc (1, imagesize);
		compare=false;
		allowCompare=false;
		repaint(compare_value);
		ReleaseMutex(heapMutex);
	}
	resetCompare();
	repaint(compare_value);
}
void grabImage()
{
	WaitForSingleObject( heapMutex, INFINITE);
	memcpy(copy_bitmap, render_bitmap_surface, imagesize);
	allowCompare=true;
	resetCompare();
	ReleaseMutex(heapMutex);
	if (lastPercent!=(int)(((double)noSimilarPixels*100)/maxPix))
		repaint(compare_value);
}
void setCompare(bool x)
{
	if (compare == x)
		return;
	if (x)
	{
		resetCompare();
		if (lastPercent!=(int)(((double)noSimilarPixels*100)/maxPix))
			repaint(compare_value);
	}
	compare = x;
}
bool getCompare()
{
	return compare;
}
int getCompareValue()
{
	if (maxPix == 0)
	{
		lastPercent=0;
		return 0;
	}
	lastPercent=(int)(((double)noSimilarPixels*100)/maxPix);
	if (lastPercent==100)
		lastPercent=lastPercent;
	return lastPercent;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											optimisation 4 - Gausian Blur
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int mask [maxsize][maxsize];
bool interpolate(int x1, int y1, int area, int Red, int Green, int Blue, int Alpha)
{
	//check if we should interpolate
	if (area <= 2 || !optimisation4 )//|| area > maxarea/10)
		return false;

	int x2 = x1+area-1;
	if (x2>=render_width)
		x2 = render_width-1;
	int y2 = y1+area-1;
	if (y2>=render_height)
		y2 = render_height-1;
	int         x ;
	int         y ;
	int         width ;
	int         height ;
	int         oldMode ;
	int         oldMode2 ;
	int         dest_width ;
	int         dest_height ;
	int         dest_xoffset ;
	int         dest_yoffset ;
	int         R ;
	int         G ;
	int         B ;
	uchar       *p ;
	uchar       *q ;
	uchar       dither ;
	HDC         hdc ;
	HDC         hdc2 ;
	RECT        rect ;
	UINT        backColor ;
	double      aspect_ratio ;
	double      screen_aspect ;
	double      preview_aspect ;

	//ExternalDisplayPlotRect (x1, y1, x2, y2, Red, Green, Blue, Alpha) ;
	if (render_priority_surface == NULL)
		return false;
	if (x1 < 0 || x1 >= render_width || x2 < 0 || x2 >= render_width)
		return false;
	if (y1 < 0 || y1 >= render_height || y2 < 0 || y2 >= render_height)
		return false;

	R = Red ;
	G = Green ;
	B = Blue ;
	for (y = y1 ; y <= y2 ; y++)
	{
	  p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x1 * 3 ;
	  for (x = x1 ; x <= x2 ; x++)
	  {
		  *p++ = Blue ;
		  *p++ = Green ;
		  *p++ = Red ;
	  }
	}
	/////////////////////////////////////////////////////////////////////////
	// INTERPOLATION    vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
	/////////////////////////////////////////////////////////////////////////
	//calculate interpolation size
	int size = area/9;
	int bwidth = size*2;
	bwidth++;
	// calculate boundries
	int xlim1 = x1-size;
	int xlim2 = x1+size;
	int xlim3 = x1+area-size;
	int xlim4 = x1+area+size;
	int ylim1 = y1-size;
	int ylim2 = y1+size;
	int ylim3 = y1+area-size;
	int ylim4 = y1+area+size;
	uchar* q2;
	//apply the mask to required cells
	//top and bottom
	for (x = xlim2 ; x < xlim3 && x>0 && x<render_width; x++)
	{
		if (ylim1>0&&ylim2<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1+1) * render_bitmap_bpl + x * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim2) * render_bitmap_bpl + x * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (y = ylim1 ; y<ylim2 && ylim1>0 ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (ylim4<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3+1) * render_bitmap_bpl + x * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim4) * render_bitmap_bpl + x * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (y = ylim3; y<ylim4 && ylim4<render_height ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}

	//left and right
	for (y = ylim2 ; y < ylim3 && y>0 && y<render_height; y++)
	{
		if (xlim1>0&&xlim2<render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim1-1) * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim2) * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (x = xlim1 ; x<xlim2 && xlim1>0 ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (xlim4<=render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim3-1) * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + xlim4 * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (x = xlim3; x<xlim4 && xlim4<render_width ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}
	//now finish off
	//top and bottom
	for (x = xlim2-1 ; xlim1<0 && x>=0 && x<render_width; x--)
	{
		if (ylim1>0&&ylim2<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1+1) * render_bitmap_bpl + x * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim2) * render_bitmap_bpl + x * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (y = ylim1 ; y<ylim2 && ylim1>0 ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (ylim4<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3+1) * render_bitmap_bpl + x * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim4) * render_bitmap_bpl + x * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (y = ylim3; y<ylim4 && ylim4<render_height ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}
	//top and bottom
	for (x = xlim3 ; xlim4>render_width && x>=0 && x<render_width; x--)
	{
		if (ylim1>0&&ylim2<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1+1) * render_bitmap_bpl + x * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim2) * render_bitmap_bpl + x * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (y = ylim1 ; y<ylim2 && ylim1>0 ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (ylim4<render_height)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3+1) * render_bitmap_bpl + x * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim4) * render_bitmap_bpl + x * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (y = ylim3; y<ylim4 && ylim4<render_height ; y++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}

	//left and right
	for (y = ylim2-1 ; ylim1<0 && y>=0 && y<render_height; y--)
	{
		if (xlim1>0&&xlim2<render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim1-1) * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim2) * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (x = xlim1 ; x<xlim2 && xlim1>0 ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (xlim4<=render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim3-1) * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + xlim4 * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (x = xlim3; x<xlim4 && xlim4<render_width ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}

	//left and right
	for (y = ylim3 ; ylim4>render_height && y>=0 && y<render_height; y++)
	{
		if (xlim1>0&&xlim2<render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim1-1) * 3 ;
			float r1 = *p++;
			float g1 = *p++;
			float b1 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim2) * 3 ;
			float r2 = *p++;
			float g2 = *p++;
			float b2 = *p++;
			float i = 1;
			for (x = xlim1 ; x<xlim2 && xlim1>0 ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r2*(i/bwidth))+(r1*((bwidth-i)/bwidth)));	
				*p++ = (int)((g2*(i/bwidth))+(g1*((bwidth-i)/bwidth)));
				*p++ = (int)((b2*(i/bwidth))+(b1*((bwidth-i)/bwidth)));
			}
		}

		if (xlim4<=render_width)
		{
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + (xlim3-1) * 3 ;
			float r3 = *p++;
			float g3 = *p++;
			float b3 = *p++;
			p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + xlim4 * 3 ;
			float r4 = *p++;
			float g4 = *p++;
			float b4 = *p++;
			float i = 1;
			for (x = xlim3; x<xlim4 && xlim4<render_width ; x++, i++)
			{
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x * 3 ;
				*p++ = (int)((r4*(i/bwidth))+(r3*((bwidth-i)/bwidth)));
				*p++ = (int)((g4*(i/bwidth))+(g3*((bwidth-i)/bwidth)));
				*p++ = (int)((b4*(i/bwidth))+(b3*((bwidth-i)/bwidth)));
			}
		}
	}


	//corners
	float ar,ag,ab,br,bg,bb,cr,cg,cb,dr,dg,db;
	for (int i = 0; i<bwidth-1; i++)
	{
		for (int j=0;j<bwidth-1;j++)
		{
			//top right
			if (xlim4<render_width && ylim1>0 && ylim1+bwidth < render_height)
			{
				//get the 4 required colours
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim3-1) * 3 ;
				ar = *p++;
				ag = *p++;
				ab = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim3+bwidth-1) * 3 ;
				br = *p++;
				bg = *p++;
				bb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1+1) * render_bitmap_bpl + (xlim3+i) * 3 ;
				cr = *p++;
				cg = *p++;
				cb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-bwidth+1) * render_bitmap_bpl + (xlim3-1+i) * 3 ;
				dr = *p++;
				dg = *p++;
				db = *p++;

				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim3+i) * 3 ;
				*p++ = (((j+1)*dr/bwidth) +((bwidth-(j+1))*cr/bwidth)+((i+1)*br/bwidth) +((bwidth-(i+1))*ar/bwidth))/2;
				*p++ = (((j+1)*dg/bwidth) +((bwidth-(j+1))*cg/bwidth)+((i+1)*bg/bwidth) +((bwidth-(i+1))*ag/bwidth))/2;
				*p++ = (((j+1)*db/bwidth) +((bwidth-(j+1))*cb/bwidth)+((i+1)*bb/bwidth) +((bwidth-(i+1))*ab/bwidth))/2;
			}

			ar=ag=ab=br=bg=bb=cr=cg=cb=dr=dg=db=255;
			//bottom right
			if (xlim4<render_width && ylim4<render_height)
			{
				//get the 4 required colours
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim3-1) * 3 ;
				ar = *p++;
				ag = *p++;
				ab = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim3+bwidth-1) * 3 ;
				br = *p++;
				bg = *p++;
				bb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3+1) * render_bitmap_bpl + (xlim3+i) * 3 ;
				cr = *p++;
				cg = *p++;
				cb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-bwidth+1) * render_bitmap_bpl + (xlim3-1+i) * 3 ;
				dr = *p++;
				dg = *p++;
				db = *p++;

				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim3+i) * 3 ;
				*p++ = (((j+1)*dr/bwidth) +((bwidth-(j+1))*cr/bwidth)+((i+1)*br/bwidth) +((bwidth-(i+1))*ar/bwidth))/2;
				*p++ = (((j+1)*dg/bwidth) +((bwidth-(j+1))*cg/bwidth)+((i+1)*bg/bwidth) +((bwidth-(i+1))*ag/bwidth))/2;
				*p++ = (((j+1)*db/bwidth) +((bwidth-(j+1))*cb/bwidth)+((i+1)*bb/bwidth) +((bwidth-(i+1))*ab/bwidth))/2;
			}

			//top left
			if (xlim1>0 && ylim1>0)
			{
				//get the 4 required colours
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim1-1) * 3 ;
				ar = *p++;
				ag = *p++;
				ab = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim1+bwidth-1) * 3 ;
				br = *p++;
				bg = *p++;
				bb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1+1) * render_bitmap_bpl + (xlim1+i) * 3 ;
				cr = *p++;
				cg = *p++;
				cb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-bwidth+1) * render_bitmap_bpl + (xlim1-1+i) * 3 ;
				dr = *p++;
				dg = *p++;
				db = *p++;

				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim1-j) * render_bitmap_bpl + (xlim1+i) * 3 ;
				*p++ = (((j+1)*dr/bwidth) +((bwidth-(j+1))*cr/bwidth)+((i+1)*br/bwidth) +((bwidth-(i+1))*ar/bwidth))/2;
				*p++ = (((j+1)*dg/bwidth) +((bwidth-(j+1))*cg/bwidth)+((i+1)*bg/bwidth) +((bwidth-(i+1))*ag/bwidth))/2;
				*p++ = (((j+1)*db/bwidth) +((bwidth-(j+1))*cb/bwidth)+((i+1)*bb/bwidth) +((bwidth-(i+1))*ab/bwidth))/2;
			}

			//bottom left
			if (xlim1>0 && ylim4<render_height)
			{
				//get the 4 required colours
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim1-1) * 3 ;
				ar = *p++;
				ag = *p++;
				ab = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim1+bwidth-1) * 3 ;
				br = *p++;
				bg = *p++;
				bb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3+1) * render_bitmap_bpl + (xlim1+i) * 3 ;
				cr = *p++;
				cg = *p++;
				cb = *p++;
				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-bwidth+1) * render_bitmap_bpl + (xlim1-1+i) * 3 ;
				dr = *p++;
				dg = *p++;
				db = *p++;

				p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - ylim3-j) * render_bitmap_bpl + (xlim1+i) * 3 ;
				*p++ = (((j+1)*dr/bwidth) +((bwidth-(j+1))*cr/bwidth)+((i+1)*br/bwidth) +((bwidth-(i+1))*ar/bwidth))/2;
				*p++ = (((j+1)*dg/bwidth) +((bwidth-(j+1))*cg/bwidth)+((i+1)*bg/bwidth) +((bwidth-(i+1))*ag/bwidth))/2;
				*p++ = (((j+1)*db/bwidth) +((bwidth-(j+1))*cb/bwidth)+((i+1)*bb/bwidth) +((bwidth-(i+1))*ab/bwidth))/2;
			}
		}
	}

	//calculate new paste area
	if (xlim1>=0)
		x1=xlim1;
	if (xlim4<render_width)
		x2=xlim4;
	if (ylim1>=0)
		y1=ylim1;
	if (ylim4<render_height)
		y2=ylim4;
	/////////////////////////////////////////////////////////////////////////
	// INTERPOLATION    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	/////////////////////////////////////////////////////////////////////////

	//only paste changes to screen if required
	if (!getScreen())
		return true;
	// WE WILL INTRODUCE CODE SO WE ALSO UPDATE PREVIEW WINDOW
	if (render_window != NULL && !IsIconic (render_window))
	{
	hdc = GetDC (render_window) ;
	oldMode = SetStretchBltMode (hdc, STRETCH_DELETESCANS) ;
	if (hPalApp)
	{
	  SelectPalette (hdc, hPalApp, false) ;
	  RealizePalette (hdc) ;
	}
	if (IsZoomed (render_window))
	{
	  //CALCULATE FOR WINDOW
	  GetClientRect  (render_window, &rect) ;
	  aspect_ratio = (double) render_width / render_height ;
	  screen_aspect = (double) rect.right / rect.bottom ;
	  if (aspect_ratio >= screen_aspect)
	  {
		dest_width = rect.right ;
		dest_height = (int) ((double) rect.right / aspect_ratio) ;
	  }
	  else
	  {
		dest_width = (int) ((double) rect.bottom * aspect_ratio) ;
		dest_height = rect.bottom ;
	  }
	  x = dest_width * x1 / render_width ;
	  y = dest_height * y1 / render_height ;
	  width = (dest_width * (x2 - x1 + 1) + render_width - 1) / render_width ;
	  height = (dest_height * (y2 - y1 + 1) + render_height - 1) / render_height ;
	  dest_xoffset = (rect.right - dest_width) / 2 ;
	  dest_yoffset = (rect.bottom - dest_height) / 2 ;
	  StretchDIBits (hdc,
					 x + dest_xoffset,
					 y + dest_yoffset,
					 width,
					 height,
					 x1,
					 render_bitmap.header.biHeight - 1 - y2,
					 x2 - x1 + 1,
					 y2 - y1 + 1,
					 render_bitmap_surface,
					 (LPBITMAPINFO) &render_bitmap,
					 DIB_RGB_COLORS,
					 SRCCOPY) ;
	}
	else
	{
	  x = x1 ;
	  y = y1 ;
	  width = x2 - x1 + 1 ;
	  height = y2 - y1 + 1 ;
	  StretchDIBits (hdc,
					 x - renderwin_xoffset,
					 y - renderwin_yoffset,
					 width,
					 height,
					 x1,
					 render_bitmap.header.biHeight - 1 - y2,
					 x2 - x1 + 1,
					 y2 - y1 + 1,
					 render_bitmap_surface,
					 (LPBITMAPINFO) &render_bitmap,
					 DIB_RGB_COLORS,
					 SRCCOPY) ;
	}
	SetStretchBltMode (hdc, oldMode) ;
	ReleaseDC (render_window, hdc) ;

	}
	//CALCULATE FOR PREVIEW
	hdc2 = GetDC (preview_window) ;
	oldMode2 = SetStretchBltMode (hdc2, STRETCH_DELETESCANS) ;
	if (hPalApp)
	{
	  SelectPalette (hdc2, hPalApp, false) ;
	  RealizePalette (hdc2) ;
	}
	GetClientRect  (preview_window, &rect) ;
	aspect_ratio = (double) render_width / render_height ;
	screen_aspect = (double) rect.right / rect.bottom ;
	if (aspect_ratio >= screen_aspect)
	{
		dest_width = rect.right ;
		dest_height = (int) ((double) rect.right / aspect_ratio) ;
	}
	else
	{
		dest_width = (int) ((double) rect.bottom * aspect_ratio) ;
		dest_height = rect.bottom ;
	}
	x = dest_width * x1 / render_width ;
	y = dest_height * y1 / render_height ;
	width = (dest_width * (x2 - x1 + 1) + render_width - 1) / render_width ;
	height = (dest_height * (y2 - y1 + 1) + render_height - 1) / render_height ;
	dest_xoffset = (rect.right - dest_width) / 2 ;
	dest_yoffset = (rect.bottom - dest_height) / 2 ;
	StretchDIBits (hdc2,
				 x + dest_xoffset,
				 y + dest_yoffset,
				 width,
				 height,
				 x1,
				 render_bitmap.header.biHeight - 1 - y2,
				 x2 - x1 + 1,
				 y2 - y1 + 1,
				 render_bitmap_surface,
				 (LPBITMAPINFO) &render_bitmap,
				 DIB_RGB_COLORS,
				 SRCCOPY) ;
	SetStretchBltMode (hdc2, oldMode2) ;
	ReleaseDC (preview_window, hdc2) ;

	return true;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											optimisation 3 - edge detection
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int unit = 8;
int edgeTolerance=unit;
void setToleranceValUp()
{
	edgeTolerance = edgeTolerance + unit;
	if (edgeTolerance>(unit*20))
		edgeTolerance=unit*20;
}
void setToleranceValDown()
{
	edgeTolerance = edgeTolerance - unit;
	if (edgeTolerance<0)
		edgeTolerance=0;
}
int getTolerance()
{
	return edgeTolerance/unit;
}
int getEdgeOffset(int x,int y, int area, int delay)
{
	//check we are within bounds
	if (x<0||y<0||x>=render_width||y>=render_height)
		return delay;
	
	//FIND CENTRAL PIXEL
	int offset;
	//calculate cx and cy to copy image pixel colour
	//make sure we work within image space
	int temp2=x+area;
	if (temp2>render_width)
	{
		offset = render_width-x;
	}
	else
		offset = area;
	int cx = x + (offset-1)/2;

	temp2=y+area;
	if (temp2>render_height)
	{
		offset = render_height-y;
	}
	else
		offset = area;
	int cy = y + (offset-1)/2;

	uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - cy) * render_bitmap_bpl + cx*3;
	int rcopy = *p++;
	int gcopy = *p++;
	int bcopy = *p++;

	int rcopy2, gcopy2, bcopy2;

	int cumulative = 0;

	int x1,x2,y1,y2;

	//top and bottom row
	y1=y-1;
	y2=y+area;
	x1=x-1;
	if (x1<0)
		x1 = 0;
	for (; x1<=x+offset && x1<render_width; x1++)
	{
		//Do top row
		if (y1>=0 && y1<render_height)
		{
			uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y1) * render_bitmap_bpl + x1*3;
			rcopy2 = *p++;
			gcopy2 = *p++;
			bcopy2 = *p++;

			if (rcopy2<rcopy)
				cumulative += (rcopy-rcopy2);
			else
				cumulative += (rcopy2-rcopy);
			if (gcopy2<gcopy)
				cumulative += (gcopy-gcopy2);
			else
				cumulative += (gcopy2-gcopy);
			if (bcopy2<bcopy)
				cumulative += (bcopy-bcopy2);
			else
				cumulative += (bcopy2-bcopy);
		}

		//Do bottom row
		if (y2>=0 && y2<render_height)
		{
			uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y2) * render_bitmap_bpl + x1*3;
			rcopy2 = *p++;
			gcopy2 = *p++;
			bcopy2 = *p++;

			if (rcopy2<rcopy)
				cumulative += (rcopy-rcopy2);
			else
				cumulative += (rcopy2-rcopy);
			if (gcopy2<gcopy)
				cumulative += (gcopy-gcopy2);
			else
				cumulative += (gcopy2-gcopy);
			if (bcopy2<bcopy)
				cumulative += (bcopy-bcopy2);
			else
				cumulative += (bcopy2-bcopy);
		}
	}

	//left and right columns
	x2=x-1;
	x2=x+area;
	y1=y;
	if (y1<0)
		y1 = 0;
	for (; y1<y+offset && y1<render_height; y1++)
	{
		//Do left column
		if (x1>=0 && x1<render_width)
		{
			uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y1) * render_bitmap_bpl + x1*3;
			rcopy2 = *p++;
			gcopy2 = *p++;
			bcopy2 = *p++;

			if (rcopy2<rcopy)
				cumulative += (rcopy-rcopy2);
			else
				cumulative += (rcopy2-rcopy);
			if (gcopy2<gcopy)
				cumulative += (gcopy-gcopy2);
			else
				cumulative += (gcopy2-gcopy);
			if (bcopy2<bcopy)
				cumulative += (bcopy-bcopy2);
			else
				cumulative += (bcopy2-bcopy);
		}

		//Do right column
		if (x2>=0 && x2<render_width)
		{
			uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y1) * render_bitmap_bpl + x2*3;
			rcopy2 = *p++;
			gcopy2 = *p++;
			bcopy2 = *p++;

			if (rcopy2<rcopy)
				cumulative += (rcopy-rcopy2);
			else
				cumulative += (rcopy2-rcopy);
			if (gcopy2<gcopy)
				cumulative += (gcopy-gcopy2);
			else
				cumulative += (gcopy2-gcopy);
			if (bcopy2<bcopy)
				cumulative += (bcopy-bcopy2);
			else
				cumulative += (bcopy2-bcopy);
		}
	}

	int temp = maxDelay - delay;
	int deficit = 0;
	//delay more if not at edge
	if (cumulative <= edgeTolerance*((area*area)+2))
	{
		deficit = ((maxDelay - delay)*2)/3;

		//int temp = area % 3;
		//if (temp == 1)
		//	deficit = (maxPix*2)/3;
		//if (temp == 2)
		//	deficit = maxPix/3;
	}

	return delay+deficit;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											optimisation 2 - randomization/delay
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double logBase(double x, double y)
{
	return log10(y)/log10(x);
}
const float constval = 10;
int timeVal = 0;
int Red=255;
int Green = 255;//testtotal;
int Blue = 255;//testtotal;
double logdenom = 1;
struct RGBcache{
	RGBcache():set(false){}
	bool set;
	int red, green, blue;
};
RGBcache colourCache [101];
int LabtoRGB(int L)
{
	if (L<0)
		L=0;
	if (L>100)
		L=100;
	if (!colourCache[L].set)
	{
		colourCache[L].set=true;

		//set the colour we are dealing with
		//blue
		/*double a = 36.879;
		double b = -51.291;*/
		//grey
		double a = 0;
		double b = 0;
		//set variables
		double varX,varY,varZ,X,Y,Z,varR,varG,varB,R,G,B;

		//first cie-lab to xyz
		varY = ((double)L+16)/116;
		varX = (a/500)+varY;
		varZ = varY - (b/200);

		if (pow (varY,3)>0.008856){
			varY = pow(varY,3);}
		else {
			varY = (varY-16/116)/7.787;}
		if (pow (varX,3)>0.008856) {
			varX = pow (varX,3);}
		else {
			varX = (varX-16/116)/7.787;}
		if (pow (varZ,3)>0.008856){
			varZ = pow (varZ,3);}
		else {
			varZ = (varZ-16/116)/7.787;}

		X = 95.047*varX;
		Y = 100*varY;
		Z = 108.88*varZ;

		//now convert to RGB
		varX = X/100;
		varY = Y/100;
		varZ = Z/100;

		varR = (varX * 3.2406) + (varY * -1.5372) + (varZ * -0.4986);
		varG = (varX * -0.9689) + (varY * 1.8758) + (varZ * 0.0415);
		varB = (varX * 0.0557) + (varY * -0.2040) + (varZ * 1.0570);

		if ( varR > 0.0031308 ) varR = 1.055 * ( pow(varR , ( 1 / 2.4 )) ) - 0.055;
		else                     varR = 12.92 * varR;
		if ( varG > 0.0031308 ) varG = 1.055 * ( pow (varG , ( 1 / 2.4 )) ) - 0.055;
		else                     varG = 12.92 * varG;
		if ( varB > 0.0031308 ) varB = 1.055 * ( pow (varB , ( 1 / 2.4 )) ) - 0.055;
		else                     varB = 12.92 * varB;

		R = varR * 255;
		G = varG * 255;
		B = varB * 255;

		//check we are not out of range
		if (R<0)
			R = 0;
		if (G<0)
			G = 0;
		if (B<0)
			B = 0;
		if (R>255)
			R=255;
		if (G>255)
			G=255;
		if (B>255)
			B=255;


		colourCache[L].red=(int)R;
		colourCache[L].green=(int)G;
		colourCache[L].blue=(int)B;
	}
	return L;
}
void priorityImage(int x1, int y1, int area, int delay)
{
	if (delay < 0)
		delay = 0;
	if (delay >= maxDelay)
		delay = maxDelay-1;

	//int diff = maxPriority-priority;
	//double top = logBase(maxPix,(double)diff);

	int val = (int)(logBase(9*maxPix,maxDelay-delay)/logdenom*100);//(int)( top/logdenom*765);

	val=100-val;

	//change contrast
	//val = 50+(50*sin (((double)val*0.031415926535897932384626433832795)-1.5707963267948966192313216916398));
	

	//check val is within range
	if (val<0)
		val = 0;
	if (val > 100)
		val = 100;

	//now that we have our value we want to use an additive colour space.
	//we will is Lab. Convertions are approximate and costly so cache values
	if (!colourCache[val].set)
		val = LabtoRGB(val);

	//sort the colours out
	Red=colourCache[val].red;
	Green=colourCache[val].green;
	Blue=colourCache[val].blue;

	//always keep alpha zero so we cant see through
	int Alpha = 0;//testtotal;
	int x2 = x1+area-1;
	if (x2>=render_width)
		x2 = render_width-1;
	int y2 = y1+area-1;
	if (y2>=render_height)
		y2 = render_height-1;
	int         x ;
	int         y ;
	int         width ;
	int         height ;
	int         oldMode ;
	int         oldMode2 ;
	int         dest_width ;
	int         dest_height ;
	int         dest_xoffset ;
	int         dest_yoffset ;
	int         R ;
	int         G ;
	int         B ;
	uchar       *p ;
	uchar       dither ;
	HDC         hdc ;
	HDC         hdc2 ;
	RECT        rect ;
	UINT        backColor ;
	double      aspect_ratio ;
	double      screen_aspect ;
	double      preview_aspect ;

	//ExternalDisplayPlotRect (x1, y1, x2, y2, Red, Green, Blue, Alpha) ;
	if (render_priority_surface == NULL)
	return ;
	if (x1 < 0 || x1 >= render_width || x2 < 0 || x2 >= render_width)
	return ;
	if (y1 < 0 || y1 >= render_height || y2 < 0 || y2 >= render_height)
	return ;
	if (x1 == x2 && y1 == y2)
	{
	//WIN_Display_Plot (x1, y1, Red, Green, Blue, Alpha,1) ;
	//return ;
	}
	R = Red ;
	G = Green ;
	B = Blue ;
	for (y = y1 ; y <= y2 ; y++)
	{
	  p = render_priority_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x1 * 3 ;
	  for (x = x1 ; x <= x2 ; x++)
	  {
		  *p++ = Blue ;
		  *p++ = Green ;
		  *p++ = Red ;
	  }
	}
	if (getScreen())
	return;
	// WE WILL INTRODUCE CODE SO WE ALSO UPDATE PREVIEW WINDOW
	if (render_window != NULL && !IsIconic (render_window))
	{
	hdc = GetDC (render_window) ;
	oldMode = SetStretchBltMode (hdc, STRETCH_DELETESCANS) ;
	if (hPalApp)
	{
	  SelectPalette (hdc, hPalApp, false) ;
	  RealizePalette (hdc) ;
	}
	if (IsZoomed (render_window))
	{
	  //CALCULATE FOR WINDOW
	  GetClientRect  (render_window, &rect) ;
	  aspect_ratio = (double) render_width / render_height ;
	  screen_aspect = (double) rect.right / rect.bottom ;
	  if (aspect_ratio >= screen_aspect)
	  {
		dest_width = rect.right ;
		dest_height = (int) ((double) rect.right / aspect_ratio) ;
	  }
	  else
	  {
		dest_width = (int) ((double) rect.bottom * aspect_ratio) ;
		dest_height = rect.bottom ;
	  }
	  x = dest_width * x1 / render_width ;
	  y = dest_height * y1 / render_height ;
	  width = (dest_width * (x2 - x1 + 1) + render_width - 1) / render_width ;
	  height = (dest_height * (y2 - y1 + 1) + render_height - 1) / render_height ;
	  dest_xoffset = (rect.right - dest_width) / 2 ;
	  dest_yoffset = (rect.bottom - dest_height) / 2 ;
	  StretchDIBits (hdc,
					 x + dest_xoffset,
					 y + dest_yoffset,
					 width,
					 height,
					 x1,
					 render_bitmap.header.biHeight - 1 - y2,
					 x2 - x1 + 1,
					 y2 - y1 + 1,
					 render_priority_surface,
					 (LPBITMAPINFO) &render_bitmap,
					 DIB_RGB_COLORS,
					 SRCCOPY) ;
	}
	else
	{
	  x = x1 ;
	  y = y1 ;
	  width = x2 - x1 + 1 ;
	  height = y2 - y1 + 1 ;
	  StretchDIBits (hdc,
					 x - renderwin_xoffset,
					 y - renderwin_yoffset,
					 width,
					 height,
					 x1,
					 render_bitmap.header.biHeight - 1 - y2,
					 x2 - x1 + 1,
					 y2 - y1 + 1,
					 render_priority_surface,
					 (LPBITMAPINFO) &render_bitmap,
					 DIB_RGB_COLORS,
					 SRCCOPY) ;
	}
	SetStretchBltMode (hdc, oldMode) ;
	ReleaseDC (render_window, hdc) ;

	}
	//CALCULATE FOR PREVIEW
	hdc2 = GetDC (preview_window) ;
	oldMode2 = SetStretchBltMode (hdc2, STRETCH_DELETESCANS) ;
	if (hPalApp)
	{
	  SelectPalette (hdc2, hPalApp, false) ;
	  RealizePalette (hdc2) ;
	}
	GetClientRect  (preview_window, &rect) ;
	aspect_ratio = (double) render_width / render_height ;
	screen_aspect = (double) rect.right / rect.bottom ;
	if (aspect_ratio >= screen_aspect)
	{
		dest_width = rect.right ;
		dest_height = (int) ((double) rect.right / aspect_ratio) ;
	}
	else
	{
		dest_width = (int) ((double) rect.bottom * aspect_ratio) ;
		dest_height = rect.bottom ;
	}
	x = dest_width * x1 / render_width ;
	y = dest_height * y1 / render_height ;
	width = (dest_width * (x2 - x1 + 1) + render_width - 1) / render_width ;
	height = (dest_height * (y2 - y1 + 1) + render_height - 1) / render_height ;
	dest_xoffset = (rect.right - dest_width) / 2 ;
	dest_yoffset = (rect.bottom - dest_height) / 2 ;
	StretchDIBits (hdc2,
				 x + dest_xoffset,
				 y + dest_yoffset,
				 width,
				 height,
				 x1,
				 render_bitmap.header.biHeight - 1 - y2,
				 x2 - x1 + 1,
				 y2 - y1 + 1,
				 render_priority_surface,
				 (LPBITMAPINFO) &render_bitmap,
				 DIB_RGB_COLORS,
				 SRCCOPY) ;
	SetStretchBltMode (hdc2, oldMode2) ;
	ReleaseDC (preview_window, hdc2) ;
}
int getDelay(int x, int y, int area)
{

	/*Depending on the current settings, will change how pixels are 
	prioritised. The delay function is factor of many things.*/

	//calculate base value. Every quality has its own unique value
	int delay = pixOffset[area-1];

	//check if edge detection is turned on
	if (optimisation3)
	{
		//this is to allow degradation of minimum value pixels. we loose unique values for pixels of max quality but this doesnt matter
		/*if (delay > maxPix)
			delay = delay - maxPix;*/
		delay = getEdgeOffset(x,y, area, delay);
	}

		
	//check if poisson is turned on
	int poisson = 0;
	if (optimisation2)
	{
		int random = rand();
		//int val = (int) log((minQueue->nItems()*constval)+1);
		//int upperbound=val*((float)1-((float)quality/maxquality));
		//upperbound++;
		//random = (random%upperbound);
		poisson += (random%maxPix);
	}
	else
	{
		poisson=(y*render_width)+x;
	}

	if (optimisation3&&delay>pixOffset[maxarea-1]);
		poisson=poisson/3;
	
	delay=delay + poisson;

	return delay;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											optimisation 1 - 'non-tree'
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//split a square pixel into 4
int splitPixel(int x, int y, int area, bool silent)
{
	//do the maths all at once to make it quicker
	int newquality=(area+2)/3;

	if (x<0||y<0||x>render_width||y>render_height)
		return newquality;

	//the split should be serialised
	WaitForSingleObject( heapMutex, INFINITE);

	//if area is <= 4 then split into individual pixels
	if (area <=4)
	{
		int temp = (area-1)/2;
		int centx = x+temp;
		int centy = y+temp;
		for (int x1 = x; x1<x+area && x1 < render_width; x1++){
			for (int y1 = y; y1<y+area && y1 < render_height; y1++)
			{
				if (silent && x1==centx && y1==centy)
					addToQueue(x1,y1,TraceType::silent,1,getDelay(x,y,1),false);
				else
					addToQueue(x1,y1,TraceType::trace,1,getDelay(x,y,1),false);
			}
		}
	}
	else
	{
		//add central pixel
		int cent = area/3;
		if (cent%2 == 1 && area%2==1 && area%3 !=0 )
			cent++;
		if (silent)
			addToQueue(x+cent,y+cent,TraceType::silent,newquality,getDelay(x,y,newquality),false);
		else
			addToQueue(x+cent,y+cent,TraceType::trace,newquality,getDelay(x,y,newquality),false);

		//some colculations
		int xright = x+area-newquality;
		int ybottom = y+area-newquality;

		//check we are not out of range
		if ( xright < render_width && ybottom < render_height)
		{
			//bottom right
			addToQueue(xright,ybottom,TraceType::trace,newquality,getDelay(x,y,newquality),false);	
		}


		//right column
		for (int ytemp = y; ytemp < ybottom && ytemp < render_height && xright<render_width; ytemp = ytemp + newquality)
		{
			addToQueue(xright,ytemp,TraceType::trace,newquality,getDelay(x,y,newquality),false);
		}

		//bottom row
		for (int xtemp = x; xtemp < xright && xtemp < render_width && ybottom<render_height; xtemp = xtemp + newquality)
		{
			addToQueue(xtemp,y+area-newquality,TraceType::trace,newquality,getDelay(x,y,newquality),false);
		}

		//top row
		for (int xtemp = x; xtemp < xright && xtemp < render_width; xtemp = xtemp + newquality)
		{
			addToQueue(xtemp,y,TraceType::trace,newquality,getDelay(x,y,newquality),false);
		}

		//left column
		for (int ytemp = y+newquality; ytemp < ybottom && ytemp<render_height; ytemp = ytemp + newquality)
		{
			addToQueue(x,ytemp,TraceType::trace,newquality,getDelay(x,y,newquality),false);
		}
	}
    ReleaseMutex(heapMutex);
        
	return newquality;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Vicon
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool getVicon()
{
	return vicon;
}

void setVicon(bool x)
{
	if (x == vicon)
		return;
	if (x)
	{
		vStartup();
		setupVicon();
	}
	else
	{
		vClose();
	}
	vicon = x;

}
int vUpdateCount = 0;
bool viconfirsttime = true;
Coordinate upperbound;
Coordinate lowerbound;
void setupVicon()
{
	viconfirsttime = true;
	upperbound.x=lowerbound.x=0;
	upperbound.y=lowerbound.y=0;
	upperbound.z=lowerbound.z=0;
}
Coordinate scale;
double tolerance=100;
void updateVicon()
{
	//count the pixels since last update
	vUpdateCount++;
	//no more than 5 frames a second
	if (/*vicon &&*/ vUpdateCount>lastPixelsPerSec/5)
	{
		Coordinate c = vFetch();
		//if this is the first meaningful fetch
		if (viconfirsttime && c.x!=0 && c.y!=0 && c.z!=0)
		{
			tolerance=100;
		}

		//check if we have moved enough
		if (c.x<lowerbound.x || c.x>upperbound.x || c.y<lowerbound.y || c.y>upperbound.y || c.z<lowerbound.z || c.z>upperbound.z )
		{
			//update bounds
			lowerbound.x = c.x-tolerance;
			lowerbound.y = c.y-tolerance;
			lowerbound.z = c.z-tolerance;
			upperbound.x = c.x+tolerance;
			upperbound.y = c.y+tolerance;
			upperbound.z = c.z+tolerance;

			////////////////////////////////
			// SORT THESE WIGHTS
			//vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv

			//transform from vicon coordinates to povray coordinates
			c.x = (c.x)/50;
			c.y = c.y/50;
			c.z = (c.z)/50;

			//apply new position
			newViewX(c.x,100);
			newViewY(c.z,100);
			newViewZ(c.y-100,100);

			setScroll(c.x,c.z,c.y-100);

			//^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
			// SORT THESE WIGHTS
			////////////////////////////////
		}

		vUpdateCount=0;
	}
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Housekeeping
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
bool setupOptimisations()
//setup dynamic 2D array that contains cell info
{
	//randomise the random number generator with system time
	srand((unsigned)time(0));

	//setup output string
	char out [100];
	wsprintf(out, "Setup picture %dx%d\n",render_width,render_height );
	OutputDebugString(out);
	//check image is not too big
	if (render_width>maxsize | render_height>maxsize)
	{
		MessageBox(NULL, "The maximum size allowed is 2500x2500\nReduce image size and try again", "Max Size", MB_ICONINFORMATION);
		return false;
	}

	//find maximum quality
    //find largest side
    int shortest = render_height;
    int longest = render_width;
    if (render_width<render_height){
		shortest = render_width;
        longest = render_height;
    }
    //we need to check that this quality is not too large
    if (longest-shortest < shortest){
            maxarea = longest;
	}else{
		maxarea = shortest;
	}
	//now calculate the largest number of pixels at each size for use in poisson function
	//calculate only once to increase efficiency
	//int val, val2;
	//for (int x = 1; x<maxquality+1; x++)
	//{
	//	val = render_width/x;
	//	if (render_width%x!=0)
	//		val++;
	//	val2 = render_height/x;
	//	if (render_width%x!=0)
	//		val2++;
	//	maxPix[x]=val*val2;
	//}
	//calculate the total number of pixels in image
	maxPix = render_width*render_height;
	//calculate offset for each quality value
	for (int x = 0; x<maxarea; x++)
	{
		pixOffset[maxarea-1-x]=x*maxPix;
	}
	//calculate maximum delay value
	maxDelay = pixOffset[0]+maxPix;
	//calculate denominator for use in priority image calculation
	logdenom = logBase(maxPix*9,(double)maxDelay);

	//set up queue
	initialiseFibHeap();
	//set tracing
	tracing = true;
	newview = true;
	refreshView = false;
	//setup 2D array and add initial pixels to queue
	if (resetFrame())
		return true;

	return false;
}
bool resetFrame()
{
	//find out if we need to setup for subdivision or not to prevent race conditions
	bool subdiv = optimisation1;
	//empty heap
	emptyQueue(true, true);
	//initiate cells
	for(int y = 0; y<render_height; y++)
	{
		for(int x = 0; x<render_width; x++)
		{
			//set to 0 if not optimisation 1
			pixelInfo [y][x].trace=TraceType::none;
			pixelInfo [y][x].area=1;
			if (!subdiv)
				addToQueue(x,y,TraceType::trace,1,getDelay(x,y,1),true);
		}
	}

	//if subdivision is off we are done
	if (!subdiv)
		return true;

	//set intial squares
	//find largest side
	shortest = render_height;
	longest = render_width;
	if (render_width<render_height){
		shortest = render_width;
		longest = render_height;
	}

	//we need to check that this quality is not too large
	if (longest-shortest < shortest)
	{
		//set the initialising pixel
		addToQueue(0,0,TraceType::trace,longest,getDelay(0,0,longest),true);
	}
	else
	{
		//the longest side is too large to fit so we will use shortest
		//calculate no. of columns
		int nowidth = render_width/shortest;

		//calculate no. of rows
		int noheight = render_height/shortest;

		char out [100];
		for (int i = 0; i <= noheight; i++)
		{
			for (int j = 0; j <= nowidth; j++)
			{
				addToQueue(i*area,j*area,TraceType::trace,shortest,getDelay(i*area,j*area,shortest),true);
			}
		}
	}
	return true;
}
void reset()
{
	optimisation1 = false;
	optimisation2 = false;
	optimisation3 = false;
	vicon = false;
	area = 1;
	repaintAll();
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Rendering functions
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//information required in drawing routine
PixelRender getPixel(int x, int y){

	updateVicon();
	viewUpdate();
	PixelRender c = getQueue();

	//display and update quality value in the options panel at most every 1 seconds
	if (lasttime<time(NULL)-1)
	{
		if(lastPixelsPerSec != pixelsPerSec)
		{
			lastPixelsPerSec = pixelsPerSec;
			repaint(average_value);
		}
		pixelsPerSec=0;
		lasttime = time(NULL);
	}

	//display quality value if it has changed
	if (c.area != area)
	{
		area = c.area;
		repaint(quality_value);
	}

	//if we are tracing or silent tracing
	if ((c.trace<=1))
	{
		pixelsPerSec++;
		//check if pixel needs subdividing
		if (c.area > 1)
		{
			if (c.x+((area-1)/2)<render_width && c.y+((area-1)/2)<render_height)
			{
				splitPixel(c.x,c.y,c.area,true);
			}
			else
			{
				splitPixel(c.x,c.y,c.area,false);
			}
		}
		else
		{
			//paint the square white as we have finished with it
			priorityImage(c.x, c.y, 1, maxDelay);
		}
	}

	return c;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Settings - Viewpoint
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//vector operations
Coordinate crossProduct(Coordinate a, Coordinate b)
{
	Coordinate c;
	c.x=a.y*b.z-a.z*b.y;
	c.y=a.z*b.x-a.x*b.z;
	c.z=a.x*b.y-a.y*b.x;
	return c;
}
double dotProduct(Coordinate a, Coordinate b)
{
	double c;
	c=(a.x*b.x) + (a.z*b.z) + (a.y*b.y);
	return c;
}
double vectorMagnitude(Coordinate a)
{
	return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}
Coordinate normaliseVec(Coordinate x, double k = 1)
{
	double magX = vectorMagnitude(x);
	if (magX == 0)
		return x;
	x.x=(x.x/magX)*k;
	x.y=(x.y/magX)*k;
	x.z=(x.z/magX)*k;
	return x;
}
Coordinate vectorAdd (Coordinate a, Coordinate b, int mult)
{
	a.x = a.x + (b.x*mult);
	a.y = a.y + (b.y*mult);
	a.z = a.z+ (b.z*mult);
	return a;
}
//other stuff
Coordinate focuspoint, workingfocuspoint;
Coordinate direction, workingdirection;
Coordinate right, sky, up;
Coordinate workingright, workingup;
double angle, workingangle;
void setSky(double x, double y, double z)
{
	sky.x=x;
	sky.y=y;
	sky.z=z;
	refreshView = true;
}
void setUp(double x, double y, double z)
{
	up.x=workingup.x=x;
	up.y=workingup.y=y;
	up.z=workingup.z=z;
	refreshView = true;
}
Coordinate getUp()
{
	Coordinate c;
	c.x=workingup.x;
	c.y=workingup.y;
	c.z=workingup.z;
	return c;
}
void setRight(double x, double y, double z)
{
	right.x=workingright.x=x;
	right.y=workingright.y=y;
	right.z=workingright.z=z;
	refreshView = true;
}
Coordinate getRight()
{
	Coordinate c;
	c.x=workingright.x;
	c.y=workingright.y;
	c.z=workingright.z;
	return c;
}
void setDirection(double x, double y, double z)
{
	direction.x = workingdirection.x = x;
	direction.y = workingdirection.y = y;
	direction.z = workingdirection.z = z;
	refreshView = true;
}
Coordinate getDirection()
{
	return workingdirection;
}
void setAngle(double x)
{
	workingangle = angle = x;
	refreshView = true;
	repaint(angle_value);
}
double getAngle()
{
	return workingangle;
}

void setFocus(double x, double y, double z)
{
	workingfocuspoint.x = focuspoint.x = x;
	workingfocuspoint.y = focuspoint.y = y;
	workingfocuspoint.z = focuspoint.z = z;
	refreshView = true;
	repaint(focus_value);
}
Coordinate getFocus()
{
	return workingfocuspoint;
}
Coordinate viewpoint, workingviewpoint,drawviewpoint;
Coordinate vecx,vecy,vecz;
double distance = 1;
int lastx,lasty,lastz;
void setView(double x, double y, double z)
{ 
	refreshView = true;
	drawviewpoint.x = workingviewpoint.x = viewpoint.x = x;
	drawviewpoint.y = workingviewpoint.y = viewpoint.y = y;
	drawviewpoint.z = workingviewpoint.z = viewpoint.z = z;

	//calculate distance between view and focus
	double distx = viewpoint.x-focuspoint.x;
	double disty = viewpoint.y-focuspoint.y;
	double distz = viewpoint.z-focuspoint.z;
	distance = sqrt(distx*distx+disty*disty+distz*distz);
	distance = distance/100;

	//calculate new viewz vector (depth)
	vecz.x=focuspoint.x-viewpoint.x;
	vecz.y=focuspoint.y-viewpoint.y;
	vecz.z=focuspoint.z-viewpoint.z;

	//calculate new the unit x vector (left right)
	//this is equal to the cross product of vecz and the unit vector in the 'up' direction
	// by default up = (0,1,0) but it can be defined differently in parse file
	//this is backwards to how me might think due to strange povray coordinate system
	Coordinate unit;
	unit.y=1;
	vecx = crossProduct(unit, vecz);

	//calculate new unit y vecor (up down)
	//we know it is perpendicular to X and Z
	//this is backwards to how me might think due to strange povray coordinate system
	// Z x X
	vecy=crossProduct(vecz, vecx);

	//now we need to make vector of our custom unit length, distance
	vecx=normaliseVec(vecx,distance);
	vecy=normaliseVec(vecy,distance);
	vecz=normaliseVec(vecz,distance);

	//setup change variables & sliders
	lastx = lasty = 0;
	lastz = -100;
	setScroll(0,0,-100);
}
Coordinate getViewpoint()
{
	return workingviewpoint;
}
Coordinate getDrawViewpoint()
{
	return drawviewpoint;
}
bool getViewChange()
{
	bool v = newview;
	newview = false;
	return v;
}
bool getViewRefresh()
{
	return refreshView;
}
void newViewX(int val, int limit)
{
	val = val*(double)limit/100;
	if (val!=lastx)
	{
		lastx = val;
		Coordinate work = vectorAdd(focuspoint,vecx,lastx);
		work = vectorAdd(work,vecy,lasty);
		drawviewpoint=vectorAdd(work,vecz,lastz);
		refreshView = true;
		repaint(view_value);
		repaint(direction_value);
		repaint(angle_value);
		repaint(focus_value);
	}
}
void newViewY(int val, int limit)
{
	val = val*(double)limit/100;
	if (val!=lasty)
	{
		lasty = val;
		Coordinate work = vectorAdd(focuspoint,vecx,lastx);
		work = vectorAdd(work,vecy,lasty);
		drawviewpoint=vectorAdd(work,vecz,lastz);
		refreshView = true;
		repaint(view_value);
		repaint(direction_value);
		repaint(angle_value);
		repaint(focus_value);
	}
}
void newViewZ(int val, int limit)
{
	val = val*(double)limit/100;
	if (val!=lastz)
	{
		lastz = val;
		Coordinate work = vectorAdd(focuspoint,vecx,lastx);
		work = vectorAdd(work,vecy,lasty);
		drawviewpoint=vectorAdd(work,vecz,lastz);
		refreshView = true;
		repaint(view_value);
		repaint(direction_value);
		repaint(angle_value);
		repaint(focus_value);
	}
}
int lastUpdateCount=0;

DWORD lastupdate = 0;
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											optimisation 5 - Plane projection
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double invalidateSize(Coordinate oldc, Coordinate newc)
{
	double distx, disty, distz, a,b,c, angle;
	double diff=1;

	//caculate distances
	distx = workingfocuspoint.x-oldc.x;
	disty = workingfocuspoint.y-oldc.y;
	distz = workingfocuspoint.z-oldc.z;
	a = sqrt(distx*distx+disty*disty+distz*distz);
	distx = workingfocuspoint.x-newc.x;
	disty = workingfocuspoint.y-newc.y;
	distz = workingfocuspoint.z-newc.z;
	b = sqrt(distx*distx+disty*disty+distz*distz);
	distx = oldc.x-newc.x;
	disty = oldc.y-newc.y;
	distz = oldc.z-newc.z;
	c = sqrt(distx*distx+disty*disty+distz*distz);

	//make positive
	if (a<0)
		a=-a;
	if (b<0)
		b=-b;
	if (c<0)
		c=-c;
	
	if (a!=0 && b!=0)
	{
		//calculate fraction of zoom
		diff = a-b;
		if (diff<0)
			diff = -diff;
		diff = (diff)/a;

		//calculate angle using cosine rule
		angle = acos((c*c - a*a - b*b)/(-a*b*2));
		diff+= (57.29578*angle)/90;

		if (diff>1)
			diff=1;
	}
	return diff;
}
struct Plane{
	Coordinate n;
	Coordinate p;
};
double sqr(double a)
{
	return a*a;
}
Coordinate addUpRight(Coordinate up, double upamount, Coordinate right, double rightamount)
{
	Coordinate a = workingfocuspoint;

	a.x = workingfocuspoint.x; 
	a.x += (up.x*upamount); 
	a.x += (right.x*rightamount);
	a.y = workingfocuspoint.y;
	a.y += (up.y*upamount);
	a.y += (right.y*rightamount);
	a.z = workingfocuspoint.z;
	a.z += (up.z*upamount);
	a.z += (right.z*rightamount);

	return a;
}
Coordinate project (Coordinate point,Plane plane, Coordinate view)
{
	Coordinate direction;
	direction.x = point.x-view.x;
	direction.y = point.y-view.y;
	direction.z = point.z-view.z;

	//calculate intersection point of line with plane
	double alpha = (dotProduct(plane.n,vectorAdd(plane.p,view,-1)))/dotProduct(plane.n,vectorAdd(point,view,-1));
	
	Coordinate project;
	//project = normaliseVec(direction,1);
	direction.x=direction.x*alpha;
	direction.y=direction.y*alpha;
	direction.z=direction.z*alpha;
	project = vectorAdd(direction, view,1);
	return project;
}
struct Matrix1D{
	Matrix1D():a(0),b(0),c(0),d(0){}
	double a,b,c,d;
};
struct Matrix2D{
	Matrix1D a,b,c,d;
};
struct Coordinate2D{
	int x,y;
};
Coordinate2D makeCoordinate2D(int x,int y)
{
	Coordinate2D c;
	c.x=x;
	c.y=y;
	return c;
}
Matrix2D getUnit()
{
	Matrix2D a;
	a.a.a=1;
	a.b.b=1;
	a.c.c=1;
	a.d.d=1;
	return a;
}
Matrix1D matrixMult(Matrix2D a,Matrix1D b)
{
	Matrix1D c;
	c.a = a.a.a*b.a+a.b.a*b.b+a.c.a*b.c+a.d.a*b.d;
	c.b = a.a.b*b.a+a.b.b*b.b+a.c.b*b.c+a.d.b*b.d;
	c.c = a.a.c*b.a+a.b.c*b.b+a.c.c*b.c+a.d.c*b.d;
	c.d = a.a.d*b.a+a.b.d*b.b+a.c.d*b.c+a.d.d*b.d;
	return c;
}
Matrix2D matrixMult(Matrix2D a,Matrix2D b)
{
	Matrix2D c;
	c.a = matrixMult(a,b.a);
	c.b = matrixMult(a,b.b);
	c.c = matrixMult(a,b.c);
	c.d = matrixMult(a,b.d);
	return c;
}
Matrix1D cToM(Coordinate c)
{
	Matrix1D a;
	a.a=c.x;
	a.b=c.y;
	a.c=c.z;
	a.d=1;
	return a;
}
Coordinate mToC(Matrix1D c)
{
	Coordinate a;
	a.x=c.a/c.d;
	a.y=c.b/c.d;
	a.z=c.c/c.d;
	return a;
}
uchar *opt5copy_bitmap;
Coordinate2D oldP1, oldP2, newP1, newP2;
void viewUpdate()
{
	//we want to count how many pixels have attempted to be drawn since last update
	lastUpdateCount++;
	//checks if the viewpoint has changed and updates fib heap if required
	if (refreshView)
	{
		double xdif = drawviewpoint.x - workingviewpoint.x;
		double ydif = drawviewpoint.y - workingviewpoint.y;
		double zdif = drawviewpoint.z - workingviewpoint.z;
		//make all differences positive
		if (xdif<0)
			xdif=-xdif;
		if (ydif<0)
			ydif=-ydif;
		if (zdif<0)
			zdif=-zdif;
		double dif = sqrt(xdif*xdif+ydif*ydif+zdif*zdif);
		if (dif>distance && (GetTickCount()-lastupdate)>=200)
		{
			lastupdate = GetTickCount();
			// we now need to reset the image to redraw
			// this should be atomic
			WaitForSingleObject( heapMutex, INFINITE);

			//make the copies required for uptimisation 5
			Coordinate oldviewpoint = workingviewpoint;
			Coordinate olddirection = workingdirection;
			Coordinate oldup = normaliseVec(workingup);
			Coordinate oldright = normaliseVec(workingright);
			Coordinate focus = workingfocuspoint;
			
			//find all new info
			updateVar();

			////////////////////////////////
			//  OPTIMISATION 5
			////////////////////////////////
			//plane projection
			if(optimisation5)
			{
				//step 2
				Plane oldplane, newplane;
				oldplane.p = newplane.p = focus;
				oldplane.n = olddirection;
				newplane.n = workingdirection;

				//step 3
				//check viewpoint is on correct side of plane
				double side = oldplane.n.x*workingviewpoint.x+oldplane.n.x*workingviewpoint.x+oldplane.n.x*workingviewpoint.x;
				if (side < 0);
				{
					//step 4
					//calculate image width
					double oldwidth = sqrt(sqr(focus.x-oldviewpoint.x)+sqr(focus.y-oldviewpoint.y)+sqr(focus.z-oldviewpoint.z));
					oldwidth = oldwidth/tan(workingangle/2);
					double newwidth = sqrt(sqr(focus.x-workingviewpoint.x)+sqr(focus.y-workingviewpoint.y)+sqr(focus.z-workingviewpoint.z));
					newwidth = newwidth/tan(workingangle/2);
					double oldheight = (oldwidth/render_width)*render_height;
					double newheight = (newwidth/render_width)*render_height;

					//step 5
					//calculate corner positions in world coordinates
					Coordinate wotl, wotr, wobl,wobr,wntl,wntr,wnbl,wnbr;
					//normalise
					Coordinate nright,nup;
					nup = (normaliseVec(workingup));
					nright = normaliseVec(workingright);
					wotl = addUpRight(oldup, oldheight, oldright, -oldwidth);
					wotr = addUpRight(oldup, oldheight, oldright, oldwidth);
					wobl = addUpRight(oldup, -oldheight, oldright, -oldwidth);
					wobr = addUpRight(oldup, -oldheight, oldright, oldwidth);
					wntl = addUpRight(nup, newheight, nright, (-newwidth));
					wntr = addUpRight(nup, newheight, nright, newwidth);
					wnbl = addUpRight(nup, -newheight, nright, (-newwidth));
					wnbr = addUpRight(nup, -newheight, nright, newwidth);

					//step 6
					//project old points into new image plane
					wotl = project (wotl, newplane, workingviewpoint);
					wotr = project (wotr, newplane, workingviewpoint);
					wobl = project (wobl, newplane, workingviewpoint);
					wobr = project (wobr, newplane, workingviewpoint);

					//step 7
					////create transformation matrix
					////7.1 subtract camera from look-at
					//Coordinate n = vectorAdd(focus, workingviewpoint,-1);
					////7.2 cross product of n and world y axis
					//Coordinate yaxis; yaxis.y=1;
					//Coordinate u = crossProduct(yaxis,n);
					////7.3 cross of u and up
					//Coordinate v = crossProduct(u,n);
					////7.4 create the matrix
					//Matrix2D transform = getUnit();
					//transform.a.a=u.x;
					//transform.a.b=u.y;
					//transform.a.c=u.z;
					//transform.a.d=-dotProduct(u,workingviewpoint);
					//transform.b.a=v.x;
					//transform.b.b=v.y;
					//transform.b.c=v.z;
					//transform.b.d=-dotProduct(v,workingviewpoint);
					//transform.c.a=n.x;
					//transform.c.b=n.y;
					//transform.c.c=n.z;
					//transform.c.d=-dotProduct(n,workingviewpoint);


					//7.1 translate to viewing coordinates
					Matrix2D t = getUnit();
					t.d.a = -workingviewpoint.x;
					t.d.b = -workingviewpoint.y;
					t.d.c = -workingviewpoint.z;

					//7.3 transform viewpoint and focus into new coordinate system
					Matrix1D f2 = cToM(focuspoint);
					//f2 = matrixMult(matrixMult(s,t),f2);
					f2 = matrixMult(t,f2);

					//7.4 rotate into yz about y
					Coordinate f = mToC(f2);
					double theta = acos(f.z/sqrt(sqr(f.x)+sqr(f.z)));
					Matrix2D r1 = getUnit();
					r1.a.a=cos(theta);
					r1.c.a=-sin(theta);
					r1.a.c=sin(theta);
					r1.c.c=cos(theta);
					//7.5 rotate about x to align with z
					f2 = cToM(f);
					Matrix1D f3 = matrixMult(r1,f2);
					f = mToC(f3);
					theta = acos(f.z/sqrt(sqr(f.y)+sqr(f.z)));
					Matrix2D r2 = getUnit();
					r2.b.b=cos(theta);
					r2.c.b=sin(theta);
					r2.b.c=-sin(theta);
					r2.c.c=cos(theta);
					//7.6 rotate into yz about y
					Matrix1D u4 = matrixMult(matrixMult(r2,r1),cToM(workingup));
					f = mToC(u4);
					theta = acos(f.y/sqrt(sqr(f.x)+sqr(f.y)));
					Matrix2D r3 = getUnit();
					r3.a.a=cos(theta);
					r3.b.a=sin(theta);
					r3.a.b=-sin(theta);
					r3.b.b=cos(theta);
					//7.7 combine to create final matrix
					Matrix2D transform = matrixMult(matrixMult(matrixMult(r3,r2),r1),t);

					//step 8 - convert all corners to viewing coordinates
					wotl = mToC(matrixMult(transform,cToM(wotl)));
					wotr = mToC(matrixMult(transform,cToM(wotr)));
					wobl = mToC(matrixMult(transform,cToM(wobl)));
					wobr = mToC(matrixMult(transform,cToM(wobr)));
					wntl = mToC(matrixMult(transform,cToM(wntl)));
					wntr = mToC(matrixMult(transform,cToM(wntr)));
					wnbl = mToC(matrixMult(transform,cToM(wnbl)));
					wnbr = mToC(matrixMult(transform,cToM(wnbr)));

					//step 9 - convert to pixel values
					Coordinate2D wntlc = makeCoordinate2D(0,0);
					Coordinate2D wntrc = makeCoordinate2D(render_width,0);
					Coordinate2D wnblc = makeCoordinate2D(0,render_height);
					Coordinate2D wnbrc = makeCoordinate2D(render_width,render_height);
					//find the scaling factor
					double scale = (render_width-1)/(wntr.x-wntl.x);
					//find pixel values of old image
					Coordinate2D wotlc = makeCoordinate2D(scale*(wotl.x-wntl.x),-scale*(wotl.y-wntl.y));
					Coordinate2D wotrc = makeCoordinate2D(scale*(wotr.x-wntl.x),-scale*(wotr.y-wntl.y));
					Coordinate2D woblc = makeCoordinate2D(scale*(wobl.x-wntl.x),-scale*(wobl.y-wntl.y));
					Coordinate2D wobrc = makeCoordinate2D(scale*(wobr.x-wntl.x),-scale*(wobr.y-wntl.y));
					
					//copy old image
					memcpy(opt5copy_bitmap, render_bitmap_surface, imagesize);

					//step 10 - the filling algorithm
					for(int y = 0; y<render_height; y++)
					{
						//check if we have anything to draw
						if ((y>=wotlc.y || y>=wotrc.y)&&(y<=woblc.y || y<=wobrc.y))
						{
							//find which lines we will intersect
							bool top = (((y>=wotlc.y && y<=wotrc.y) || (y<=wotlc.y && y>=wotrc.y))&&(wotlc.y!=wotrc.y));
							bool left = (y>=wotlc.y && y<=woblc.y);
							bool bottom = (((y>woblc.y && y<wobrc.y) || (y<wobrc.y && y>woblc.y))&&(woblc.y!=wobrc.y));
							bool right = (y>=wotrc.y && y<wobrc.y);

							//calculate how far along line we are and find pairs of points
							double d;
							bool set=false;
							if (left)
							{
								set = true;
								//find distance along line
								double denom = wotlc.y-woblc.y;
								d = (wotlc.y - y) / denom;
								//find x pos in new coords
								denom = woblc.x-wotlc.x;
								newP1.x=wotlc.x+(d*denom);
								newP1.y=y;
								//find position in oldcoords
								oldP1.x=0;
								oldP1.y=d*(render_height-1);
							}

							if (top)
							{
								set = true;
								//find distance along line
								double denom = wotrc.y-wotlc.y;
								if (denom<0)
									denom = -denom;
								d = (y - wotlc.y) / denom;
								//find x pos in new coords
								denom = wotrc.x-wotlc.x;
								if(set)
								{
									newP2.x=wotlc.x+(d*denom);
									newP2.y=y;
									//find position in oldcoords
									oldP2.x=d*(render_width-1);
									oldP2.y=0;
								}
								else
								{
									newP1.x=wotlc.x+(d*denom);
									newP1.y=y;
									//find position in oldcoords
									oldP1.x=d*(render_width-1);
									oldP1.y=0;
								}
							}

							if (bottom)
							{
								//find distance along line
								double denom = wobrc.y-woblc.y;
								d = (y-woblc.y) / denom;
								//find x pos in new coords
								denom = wobrc.x-woblc.x;
								if (denom<0)
									denom=-denom;
								if(set)
								{
									newP2.x=woblc.x+(d*denom);
									newP2.y=y;
									//find position in oldcoords
									oldP2.x=d*(render_width-1);
									oldP2.y=(render_height-1);
								}
								else
								{
									newP1.x=woblc.x+(d*denom);
									newP1.y=y;
									//find position in oldcoords
									oldP1.x=d*(render_width-1);
									oldP1.y=(render_height-1);
								}
							}

							if (right)
							{
								//find distance along line
								double denom = wotrc.y-wobrc.y;
								d = (wotrc.y - y) / denom;
								//find x pos in new coords
								denom = wobrc.x-wotrc.x;
								newP2.x=wotrc.x+(d*denom);
								newP2.y=y;
								//find position in oldcoords
								oldP2.x=(render_width-1);
								oldP2.y=d*(render_height-1);
							}

							//order points
							if (newP1.x>newP2.x)
							{
								Coordinate2D temp = newP2;
								newP2 = newP1;
								newP1 = temp;

								temp  = oldP2;
								oldP2 = oldP1;
								oldP1 = temp;
							}
							//do calculation just once
							int xdif = oldP2.x-oldP1.x;
							int ydif = oldP2.y-oldP1.y;
							int denom = newP2.x-newP1.x;
							for(int x = 0; x<render_width; x++)
							{
								if (x>=newP1.x && x<=newP2.x && denom!=0)
								{
									double frac = ((double)x-newP1.x)/denom;
									int oldx = (frac*xdif)+oldP1.x;
									int oldy = (frac*ydif)+oldP1.y;
									uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x*3;
									uchar* p2 = opt5copy_bitmap + (render_bitmap.header.biHeight - 1 - oldy) * render_bitmap_bpl + oldx*3;
									*p++ = *p2++;
									*p++ = *p2++;
									*p++ = *p2++;
								}
								else
								{
									uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x*3;
									uchar* p2 = opt5copy_bitmap + (render_bitmap.header.biHeight - 1) * render_bitmap_bpl;
									*p++ = *p2++;
									*p++ = *p2++;
									*p++ = *p2++;
								}
							}
						}
						else
						{
							for(int x = 0; x<render_width; x++)
							{
								uchar* p = render_bitmap_surface + (render_bitmap.header.biHeight - 1 - y) * render_bitmap_bpl + x*3;
								uchar* p2 = opt5copy_bitmap + (render_bitmap.header.biHeight - 1) * render_bitmap_bpl;
								*p++ = *p2++;
								*p++ = *p2++;
								*p++ = *p2++;
							}
						}
					}
				}
				resetCompare();
				repaint(preview_window);
				InvalidateRect (render_window, NULL, TRUE);
			}
			resetCompare();
			repaint(preview_window);
			//find out if we need to setup for subdivision or not to prevent race conditions
			if(!optimisation1)
			{
				//initiate cells
				for(int y = 0; y<render_height; y++)
				{
					for(int x = 0; x<render_width; x++)
					{
						if (pixelInfo [y][x].trace<=TraceType::silent)
						{
							pixelInfo [y][x].trace = TraceType::trace;
							pixelInfo [y][x].area = 1;
						}
						else
						{
							addToQueue(x,y,TraceType::trace,1,getDelay(x,y,1),false);
						}
					}
				}
			}
			else
			{
				//empty heap
				emptyQueue(true, false);

				//calculate size of subdivision
				double frac = invalidateSize(oldviewpoint, workingviewpoint)*maxarea;
				int interval = (int) frac;
				// if interval is outside of bounds, correct
				if (interval<1)
					interval = 1;
				if (interval>maxarea)
					interval = maxarea;
				//we should only invalidate at this quality if greater than current render level
				if (interval<area)
					interval = area;
				//setup ints
				int yinterval = 0;
				int xinterval = 0;
				//invalidate cells
				for(int y = 0; y<render_height; y++)
				{
					xinterval = 0;
					for(int x = 0; x<render_width; x++)
					{
						//set to 0 if not optimisation 1
						pixelInfo [y][x].trace=TraceType::none;
						pixelInfo [y][x].area=1;
						if(x == xinterval & y == yinterval)
						{
							addToQueue(x,y,TraceType::trace,interval,getDelay(x,y,interval),false);
						}
						if(x > xinterval)
						{
							xinterval += interval;
						}
						if(y > yinterval)
						{
							yinterval += interval;
						}
						char out [100];
						wsprintf(out, "interval %d\n", interval);
						if (false)
						{
							OutputDebugString(out);
						}
					}
				}
			}

			//reset these variables
			lastUpdateCount = 0;
			refreshView=false;
			newview=true;
			repaint(view_value);
			repaint(direction_value);
			repaint(angle_value);
			repaint(focus_value);
			ReleaseMutex(heapMutex);
		}
	}
}
void updateVar()
{
	workingviewpoint = drawviewpoint;
	//we need to update all the other variables too
	//first calculate some lengths
	double dlength = vectorMagnitude(direction);
	double ulength = vectorMagnitude(up);
	double rlength = vectorMagnitude(right);
	Coordinate tempv = crossProduct(up, direction);
	double Handedness = dotProduct(tempv, right);

	//direction = focus-viewpoint
	workingdirection.x = workingfocuspoint.x - workingviewpoint.x;
	workingdirection.y = workingfocuspoint.y - workingviewpoint.y;
	workingdirection.z = workingfocuspoint.z - workingviewpoint.z;

	// Check for zero length direction vector.
	if (workingdirection.x+workingdirection.y+workingdirection.z < 1.0e-10 && workingdirection.x+workingdirection.y+workingdirection.z > -1.0e-10)
	{
		workingdirection.x = workingdirection.y = workingdirection.z = 1;
	}

	workingdirection = normaliseVec(workingdirection,1);

	//we also need to update 'up' and 'right' vectors

	workingright = crossProduct(sky, workingdirection);

	workingright = normaliseVec(workingright, 1);
	workingup = crossProduct(workingdirection,workingright);
	workingdirection = normaliseVec(workingdirection, dlength);

	if (Handedness > 0.0)
	{
		workingright = normaliseVec (workingright, rlength);
	}
	else
	{
		workingright = normaliseVec (workingright, -rlength);
	}

	workingup = normaliseVec(workingup, ulength);


	repaint(view_value);
	repaint(direction_value);
	repaint(angle_value);
	repaint(focus_value);
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Settings - Quality
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int getArea()
{
	return area;
}
int getAvgTime()
{
	return lastPixelsPerSec;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//											Settings - Optimisations
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
bool getOptimsation1()
{
	return optimisation1;
}
bool getOptimsation2()
{
	return optimisation2;
}
bool getOptimsation3()
{
	return optimisation3;
}
bool getOptimsation4()
{
	return optimisation4;
}
bool getOptimsation5()
{
	return optimisation5;
}
void setOptimsation1(bool x)
{
	if (optimisation1 == x)
		return;
		//set bool
	optimisation1 = x;
	if (!x)
	{
		setOptimsation2(false);
		setOptimsation4(false);
	}
	if (!tracing)
		return;
	working = true;
	repaint(workinglabel);
	// the whole operation needs to be atomic
	WaitForSingleObject( heapMutex, INFINITE);
	if (x)
	{
		//switch subdivision on
		//empty queue but do not empty table
		emptyQueue(false, false);
		//go through every element in image looking for undrawn
		for(int y = 0; y<render_height; y++)
		{
			for(int x = 0; x<render_width; x++)
			{
				if (pixelInfo [y][x].trace<=1)
				{
					//if undrawn, find largest area box possible
					int area = 0;
					bool error = false;
					while (!error)
					{
						area++;
						//check we have not gone off edge of image
						if (x+area>render_width || y+area>render_height)
							error = true;
						//check column
						for (int y2 = y; y2<y+area+1&&!error;y2++)
						{
							if (pixelInfo[y2][x+area].trace==TraceType::none)
								error = true;
						}
						//check row
						for (int x2 = x; x2<x+area+1&&!error;x2++)
						{
							if (pixelInfo[y+area][x2].trace==TraceType::none)
								error = true;
						}
					}
					// now update area to represent new layout
					// first unset all pixels in area
					for (int y2 = y; y2<y+area;y2++)
					{
						for (int x2 = x; x2<(x+area);x2++)
						{
							pixelInfo [y2][x2].trace=TraceType::none;
							pixelInfo [y2][x2].area=1;
							pixelInfo [y2][x2].delay=0;
						}
					}
					//finaly set pixel and add to heap
					int val = getDelay(x,y,area);
					minQueue->insert((y*render_width) + x,val);
					pixelInfo [y][x].trace=TraceType::trace;
					pixelInfo [y][x].area=area;
					pixelInfo [y][x].delay=val;
					priorityImage(x,y,pixelInfo [y][x].area,val);
				}
			}
		}
	}
	else
	{
		//switch subdivision off by subdividing all pixels to 1
		backupInfo();
		//empty queue but do not empty table
		emptyQueue(false, false);
		for(int y = 0; y<render_height; y++)
		{
			for(int x = 0; x<render_width; x++)
			{
				//set to 0 if not optimisation 1
				if (tempInfo [y][x].trace<=1)
				{
					for(int y2 = y; y2<(y+tempInfo [y][x].area); y2++)
					{
						for(int x2 = x; x2<(x+tempInfo [y][x].area); x2++)
						{
							pixelInfo [y2][x2].delay = getDelay(x2,y2,1);
							minQueue->insert((y2*render_width) + x2,pixelInfo [y2][x2].delay);
							pixelInfo [y2][x2].trace=TraceType::trace;
							pixelInfo [y2][x2].area=1;
							priorityImage(x2,y2,pixelInfo [y2][x2].area,pixelInfo [y2][x2].delay);
						}
					}
					pixelInfo [y][x].trace = tempInfo [y][x].trace;
				}
			}
		}
	}
	working = false;
	repaint(workinglabel);
	ReleaseMutex(heapMutex);
}
void setOptimsation2(bool x)
{
	if (optimisation2 == x)
		return;
	optimisation2 = x;
	if (x)
		setOptimsation1(true);
	if (!tracing)
		return;
	working = true;
	repaint(workinglabel);
	// the whole operation needs to be atomic
	WaitForSingleObject( heapMutex, INFINITE);
	//empty queue but do not empty table
	emptyQueue(false, false);
	for(int y = 0; y<render_height; y++)
	{
		for(int x = 0; x<render_width; x++)
		{
			//set to 0 if not optimisation 1
			if (pixelInfo [y][x].trace<=1)
			{
				int delayval = getDelay(x,y,pixelInfo [y][x].area);
				minQueue->insert((y*render_width) + x,delayval);
				priorityImage(x,y,pixelInfo [y][x].area,delayval);
			}
		}
	}
	ReleaseMutex(heapMutex);
	working = false;
	repaint(workinglabel);
}
void setOptimsation3(bool x)
{
	if (optimisation3 == x)
		return;
	optimisation3 = x;
	if (!tracing)
		return;
	working = true;
	repaint(workinglabel);
	// the whole operation needs to be atomic
	WaitForSingleObject( heapMutex, INFINITE);
	//empty queue but do not empty table
	emptyQueue(false, false);
	for(int y = 0; y<render_height; y++)
	{
		for(int x = 0; x<render_width; x++)
		{
			//set to 0 if not optimisation 1
			if (pixelInfo [y][x].trace<=1)
			{
				int delayval = getDelay(x,y,pixelInfo [y][x].area);
				minQueue->insert((y*render_width) + x,delayval);
				priorityImage(x,y,pixelInfo [y][x].area,delayval);
			}
		}
	}
	ReleaseMutex(heapMutex);
	working = false;
	repaint(workinglabel);
}
void setOptimsation4(bool x)
{
	if (optimisation4 == x)
		return;
	optimisation4 = x;
	if (x)
		setOptimsation1(true);
	if (!tracing)
		return;
	working = true;
	repaint(workinglabel);
	// the whole operation needs to be atomic
	working = false;
	repaint(workinglabel);
}
void setOptimsation5(bool x)
{
	if (optimisation5 == x)
		return;
	optimisation5 = x;
}
bool getTrace()
{
	return tracing;
}